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Abstract: In this paper a forward solver software for the time domain and 
the CW domain based on the Born approximation for simulating the effect 
of small localized fluorophores embedded in a non-fluorescent biological 
tissue is proposed. The fluorescence emission is treated with a mathematical 
model that describes the migration of photons from the source to the 
fluorophore and of emitted fluorescent photons from the fluorophore to the 
detector for all those geometries for which Green's functions are available. 
Subroutines written in FORTRAN that can be used for calculating the 
fluorescent signal for the infinite medium and for the slab are provided 
with a linked file. With these subroutines, quantities such as reflectance, 
transmittance, and fluence rate can be calculated. 
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1. Introduction 

Fluorescence spectroscopy and imaging using near infrared light are methods of analysis that 
have become emerging and promising techniques when applied in medicine [1-7]. Applica- 
tions rely on measurements of the steady-state emission spectrum when continuous wave (CW) 
sources are used, on measurements of the time-resolved (TR) fluorescence signal when pulsed 
laser sources are employed (time domain), or on measurements of frequency domain signals 
(amplitudes and phase-shift measurements) when modulated laser sources are employed [8]. 
The development of fluorescence applications has been accompanied by the elaboration of 
several forward solvers in order to obtain the fluorescence signal. In the literature we have 
examples of forward solvers based on solutions of the diffusion equation (DE) [9-12], on the 
random walk theory [13], on the spherical harmonics approximation [14], and on numerical so- 
lutions of the radiative transfer equation (RTE) obtained using the Monte Carlo method [15, 16] 
or with the finite-difference discrete-ordinates method [17]. The Born approximation applied 
to the DE has also been used for measurements in the frequency domain [18, 19]. For the in- 
verse problem O'Leary et al. [18] used a first-order perturbation model for the diffuse photon 
density waves emitted from a localized fluorophore inside an infinite medium. Ntziachristos 
and Weissleder [19] used a normalized Born approximation of the DE to develop an algorithm 
for the tomographic reconstruction of small fluorescent targets inside an infinite medium. Li 
et al. [20] have presented solutions for the fluorescent diffuse photon density waves of a flu- 
orescent target embedded inside an infinite or a semi-infinite fluorescent diffusive medium. 
These solutions for the frequency domain are not analytical in a proper sense since numerical 
evaluations of integrals are required. 

We propose a perturbative time domain and the CW domain forward solver software based 
on the Born approximation in order to simulate the effect of localized fluorophores embedded 
in a non fluorescent biological tissue. The fluorescence emission is treated with a mathematical 
model that is capable of describing the migration of photons from the source to the fluorophore 
and of emitted fluorescent photons from the fluorophore to the detector, for all those geometries 
for which Green's functions can be calculated. The method is implemented here for the infinite 
medium and for the slab for which some subroutines written in FORTRAN are provided with 
the linked file (Media 1). Quantities such as reflectance, transmittance, and fluence rate can 
be calculated by using these subroutines and the software package distributed with a recently 
published book [21]. The forward solver software has been developed for the time domain and 
the CW domain, however guidelines that can extend its use also to the frequency domain have 
been provided as well. 

The software presented in this paper is a tool for photon migration studies that can be utilised 
when a fluorescent inclusion is present in tissue. For instance, this is the case of optical molec- 
ular imaging that uses near-infrared fluorescent probes [14, 17]. In optical molecular imaging, 
a fluorescent biochemical marker is injected into tissue and will emit near-infrared light upon 
excitation from an external light source. By using measurements of the light intensity on the 
tissue surface, it is possible to determine the spatial concentration distribution of the marker 
inside tissue by recording tomographical data sets and employing appropriate tomographic im- 
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age reconstruction schemes [22]. The inverse fluorescent source problem in optical molecular 
imaging has been solved by using the DE [19, 23, 24] or the RTE [ 17, 25] as a forward model. 
Thus, the availability of fast and precise forward solvers is fundamental for the inverse problem 
in optical molecular imaging. 

The theory is presented in section 2; the software is described in section 3 and in the body 
text of the subroutines provided with the linked file (Media 1). In section 3, some examples 
of use of the software are also provided together with a comparison with the results of Monte 
Carlo simulations. Lastly, in section 4 a summary of the paper and a discussion of possible 
further developments of the software are provided. 

2. Theory 

The forward solver which we propose considers the effects of fluorescence on photon migration 
by the probability of having fluorescent photons emitted inside the medium. The fluorescent 




Fig. 1 . Diagram of a photon path from the source to the fluorophore and from the fluo- 
rophore to the detector. 

signal is affected by multiple scattering at both the excitation and the emission wavelengths. In 
Fig. 1, a diagram shows that the excitation light at wavelength X is an isotropic point source 
at r s , and that the fluorescent light at wavelength X e , which is emitted from a fluorophore at 
r , is collected at point r. The absorption coefficient, the reduced scattering coefficient, and the 
refractive index of the background medium at X are denoted, respectively, as jj. a , jj. s , and «, and 
at X e as ji ae , jj, se , and n e . We assume that a fluorophore is distributed inside a volume V' with 
absorption coefficient and reduced scattering coefficient jj, a f and jj, s ^ at X, and \l a f e and \l s f e at 
X e . 

A similarity can be observed between the behavior of an absorbing inclusion and that of 
a fluorescent inclusion when these are inside a scattering medium. Both of them subtract en- 
ergy from propagation in accordance to the phenomena of photons absorption. Therefore, the 
expression of the perturbed fluence due to an absorbing defect obtained with the Born approx- 
imation (see, for instance, chapter 7 in [21]) can be used to describe the photons absorbed by 
a fluorescent inclusion. The signal at emission wavelength can then be calculated by using the 
probability function that the fluorophore re-emits the absorbed photons. 

2.1. Perturbative model based on the DE 

Light propagation through biological tissue can be described with the time-dependent DE [21]. 
Therefore, the phenomena of fluorescence can be described by coupled diffusion equations at 
excitation and emission wavelengths that are respectively, [10]: 

•fia+fiaf-DV 2 



v dt 



*(r,/)=«(r,/), (1) 
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1 d 

Ve dt 



<S> e (r,t)=q e (r,t), 



(2) 



where D = 1 / [3 (jJ, s + Mi/)] an d D e = l/[3 (ji se + fi se f)] are the diffusion coefficients at excitation 
and emission wavelengths. Equations (1) and (2) are coupled by the excitation fluence, <t>, which 
determines the source term q e (r,t) of Eq. (2). In particular, the probability dPf that fluorescent 
photons are generated in the time interval dt' at time t' , in the wavelength interval dX e around 
X e , and in the volume element af 3 r' at r , is 

dP f (r , t) = G(r s , r',/i a + \l a f, \l s + «, t Va/ 1 ? , X e )d 3 r'dX e dt ' , (3) 

where r](X,X e ) is the quantum yield or quantum efficiency per unit wavelength for emission 
at X e given excitation at X, and G is the Green's function of the medium for the fluence at the 
excitation wavelength. Therefore, the coupling of Eqs. (1) and (2) is given by the term 



q e (r',t) = —7](X,X e )dX e [ exp [- - — — 1 G(r s ,r',ju„ +Haf,fl s + n' sf ,n,t')dt, 
X JO L X J J 



(4) 



where X is the fluorescent lifetime of the fluorophore (i.e. the average time that the fluorophore 
spends in the excited state prior to returning to the ground state). By substituting Eq. (4) in Eq. 
(2), the following equation at emission wavelength is obtained: 



\_d_ 

V e d7 



-yLae + H-afe-D e V 2 



(t -t 



3> e (r,/)-^T^ e / exp -i G(r s ,r ,...,t )dt =0, 

X Jo L X i 

(5) 

with r]^ e = r](X,X e )dX e quantum efficiency of the fluorophore in the wavelength range dX e 
around X e , i.e. the ratio of the number of photons emitted in the wavelength range dX e to the 
number of absorbed photons. According to Green's function theorem, the fluence rate at emis- 
sion wavelength, <t> e (r,f ), due to the fluorescent photons generated in the interval dX e , is 



4> e (r , t ) = Jo J y , q e (r ,t')G e (r',r, \i ae + li afe , + n' sfe ,n e ,t -t')dr dt 



(6) 



with G e Green's function of the medium for the fluence at the emission wavelength. We stress 
that the Green's functions G and G e are the Green's function of the same medium, so that the 
different notation is maintained only to distinguish between the excitation and the emission 
wavelengths. If all the fluorescence photons were promptly emitted, the fluence at emission 
wavelength would simply be given by 

<t> eQ (r , t ) = Ha/ri^e Jo I v ' G(r s , r , t)G e (r , r, ...,t - 1 )dx dt . (7) 
Equation (6) is thus the convolution of <i> e o with the exponential decay of the fluorophore 



<5 e (r,f) =®eo(v,t)*-exp(-t/*)- 

X 



(8) 



It can be noted that the integral of Eq. (7) is similar to the expression of the perturbed fluence 
due to an absorbing defect placed inside a turbid medium obtained using the Born approxi- 
mation (see, for instance, chapter 7 in [21], Eq. (7.11) for the DE or Eq. (7.63) for the RTE). 
Therefore, the Born approximation can also be used to calculate the fluorescence signal. That 
is, if the fluorophore occupies a small volume, and if at excitation wavelength it perturbs the 
fluence rate of the background medium of a small quantity then the following approximations 
can be made: 



G(r s ,r ,n a + Li af ,Li s + Li sf ,n,t ) ~ G(r s ,r ,^ a ,^ s ,n,t ) 

G e (r , r, [i a< , + ii aef , n' se + fj,' sef , n e , t - 1 ' ) ~ G e (r , r , fi ae , n' se , n e , t - 1 ) 

G(r s ,r ,n a ,^' s ,n,t) ^G(r s ,r" ,^ a ,n' s ,n,t), Vr',r" eV 

G e (r ,r,il ae ,il se ,n e ,t -t) ~ G e {r ,r,ii aei jx se ,n e ,t -t), Vr',r" e V 



(9) 
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and the calculation of the TR expressions for the fluence <i> e o (r,f) and <P e (r,f) can simply be 
obtained as 



*eo(r,t) = Vafri^eV' foG(r s ,r\p a ,^,n,t')G e (r ,r,Liae,Li' se ,n e ,t -t')dt , (10) 



xG f (r ,r,jj, ae ,jj, se ,n e ,t — f jar ar . 



(11) 



These equations are perturbative solutions for a fluorescence emission in which the integral on 
V is omitted, and are valid only for small fluorescent inclusions. 

Equations (10) and (11) were obtained using the DE, but can also be obtained by using 
the RTE [21]. We stress that they can be calculated for any geometry for which analytical or 
numerical Green's functions for the DE or the RTE are available. For G and G e the solutions and 
the software presented in [21] can be used and Eqs. (10) and (11) can be therefore calculated 
for all the geometries analyzed in [21]. 

For the sake of clarity we stress that, for the theory presented in this section as well as the 
models proposed in [18, 19], the Born approximation is applied twice: i.e., a first time at excita- 
tion wavelength, and a second time at emission wavelength. In fact, both the Green's functions 
G and G e are approximated to that of the background medium. Therefore, these models are 
similar but not exactly equal to the Born approximation. 

The following model is expected to fail when the fluorophore significantly perturbs the flu- 
ence at excitation wavelength and thus when the approximations of Eqs. (9) are no longer valid 
and so the actual expressions of G and G e should be used. In this case, a Green's function of the 
medium obtained with the high order perturbation theory could be used to improve the forward 
solver [26-28]. 

The theory has been developed for time-dependent sources. It is easy to extend this approach 
to CW sources. By eliminating the time dependence in Eq. (1), a CW perturbative solution 
similar to Eq. (10) can be written as follows: 

3> e (r) =^ af T]^ e v'G(r s ,r',^ a ,^ s ,n)G e (r',r,^ ae ,^ se ,n e ). (12) 

Last, we considered the frequency-domain: When modulated sources with a DC and a AC 
component were considered [29], the source term of Eq. (1) could be written as 

q(r,t) =qoDc(r)+qoAc{r)esp(-iatt), (13) 

where CO = 2k f is the pulsation and / is the frequency of modulation, qooc an d tfOAC are the 
DC and AC parts of the source strength. The solution for the fluence can be written as the sum 
of the DC and AC solutions, i.e., <£>(r,/) = <J>Dc( r ) +®Ac( r , t )■ The AC solution of the DE with 
the source term given by Eq. (13) must be oscillating at the same pulsation CO as the source. 
Thus, the AC solution of Eq. (1) will have the form: 

4> AC (r,r ) = *Acoexp(-z(Bf). (14) 

If we substitute the fluence in Eq. (1) with Eq. (14), we obtain that <t>Aco can be described by 
the equation 

[-{ico)/v + H a + ll a f-DV 2 ] 4> AC0 (r) =q 0 (r). (15) 

Equation (15) does not have a time-dependent term, and it is a CW diffusion equation. Thus, 
the general solution for <&aco can t> e obtained by using the CW solution, provided that the 
coefficient fi a is replaced with the term [ju fl — {ico)/v\. 
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2.2. Hybrid approach for the slab 

The theory described above is obtained by using the DE and is not restricted to any particu- 
lar geometry. We should mention that Eqs. (11) and (12) can also be obtained by using the 
RTE [21]. In general, solutions based on the DE are affected by the intrinsic approximations 
of this theory, which show different consequences for steady-state and time-dependent sources. 
For steady-state sources, solutions based on the DE provide an insufficient description of pho- 
tons detected at short distances [21]. For time-dependent sources, DE solutions provide an 
insufficient description of early received photons [21]. The limitations of the DE are a seri- 
ous problem in modeling photon migration, especially when small volumes of diffusive media 
are considered. To overcome these limitations, solutions that are an improvement of the DE 
solutions [30-34], hybrid solutions [35-37], or RTE solutions [38-40] have been provided. 

To evaluate the fluorescence signal with Eqs. (11) and (12), analytical expressions for the TR 
and CW Green's functions at excitation and emission wavelengths are required. The Green's 
function for the slab is calculated by using the hybrid approach described in [21], which makes 
use of RTE solutions for the infinite medium, together with the DE framework approach, which 
is based on the method of images. For the time domain (Eq. (11)), the solutions for the infinite 
medium that can be used are the solution of the DE, the solution of the telegrapher equation, 
and the solution of the RTE for isotropic scattering proposed by Paasschens [38]. For the CW 
domain (Eq. (12)), these solutions for the infinite medium are considered: the solution of the 
DE and the improved solution of the DE (IDE) for isotropic scattering proposed by Graff and 
Rinzema [32]. The solution by Graff and Rinzema will be denoted inside the software section 
as IDE solution in the CW domain. 

The limitations of the solutions implemented for the slab, i.e. of Eqs. (11) and (12) calcu- 
lated using the hybrid approach adopted in [21], can be ascribed to: I) approximations due to 
the assumptions contained in Eqs. (9) that are similar to the Born approximation; II) approx- 
imations of the hybrid Green's functions for the slab described in [21]. Improvements on this 
theory could be obtained by implementing: A) a better approximation for the Green's function 
describing the fluence at the fluorophore site, which could be obtained by using high-order so- 
lutions of the perturbation theory [26,27,42] or exact RTE solutions [39]; B) the whole integral 
on the volume V . 

3. Software and examples of use 

With this paper (see Media 1), we provide some subroutines written in FORTRAN 
that together with the subroutines described in [21] can be used for calculating Eqs. 
(11) and (12) for the infinite medium and for the homogeneous slab. The Homoge- 
neous_Fluence_Fluorescence Infinite subroutine is the basis for calculating the fluorescence 
in the infinite medium. The calculation can be implemented in the time-domain by using 
Eq. (11) and in the CW domain by using Eq. (12). The subroutines necessary for the cal- 
culation are RTE -Fluence Infinite -TR (solution by Paasschens [38]), DE -Fluence Infinite -TR, 
DE_Fluence Infinite _CW, and IDE _Fluence Infinite _CW (solution by Graaff and Rinzema 
[32]). The code of these subroutines can be found in [21] or in the linked file (Media 1). The Ho- 
mogeneous _Fluence -Fluorescence Slab subroutine is the basis for calculating the fluorescence 
in the slab. The subroutines necessary for the calculation are the same mentioned above. 

Some instructions for the use of the Homogeneous -Fluence -Fluorescence Infinite and Ho- 
mogeneous _Fluence -Fluorescence Slab subroutines are provided within the body of the text 
on the subroutines. The reflectance, R(p,t), that is due to the fluorophore can be calculated by 
using Fick's law or the partial current boundary condition [21]: 

fi(p,f)=Z)|-4>(p,z = 0,f), (16) 
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R(p,t) = ^*(p,z = 0,t), (17) 

with D being the diffusion coefficient, A the coefficient for the Fresnel Reflections [21], and <t> 
the fluorescent fluence calculated with the proposed subroutines. 

We stress that when Eqs. (11) and (12) are calculated using RTE Green's functions for the 
infinite medium, the limit of the calculated solution tends to become exact as the volume of the 
fluorescent inclusion approaches to zero. Therefore, the proposed forward solver can be exact 
for small fluorescent inclusions. By making use of the recent RTE Green's function obtained by 
Liemert and Kienle [40] or by Machida et al. [39], it is possible to build an exact forward solver 
for small localized fluorescent inclusions even for anisotropic scattering. Last, we stress that 
subroutines for calculating the fluorescence signal emitted by a localized fluorescent inclusion 
placed inside a two-layer medium, a sphere, a parallelepiped, and a cylinder, can be obtained 
by using the software package provided in [21] and by implementing for these geometries a 
similar approach to that used for the infinite medium and the slab. 
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Fig. 2. Fluorescence TR reflectance signal received from a slab 2000 mm thick at p =15 
mm. The absorption and the reduced scattering coefficient of the background medium are 
equal to 0.01 and 1. mm -1 , respectively, at both the emission and the excitation wave- 
lengths. The fluorophore is a sphere centered at (7.5, 0, 7.5) mm and with a volume V' = 
1.25 mm 3 . The absorption coefficient of the fluorophore at excitation wavelength is \L a f = 
0.002 mm~'. The refractive index of the medium and of the external is 1.4. 
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In order to cite an example of use of the proposed software, we consider a fluorophore em- 
bedded in a slab 2000 mm thick at p =15 mm. The fluorophore is assumed by a sphere centered 
at (7.5, 0, 7.5) mm with a volume V 1 = 1.25 mm 3 . The absorption and the reduced scattering 
coefficient of the background medium are equal to 0.01 and 1. mm -1 , respectively, at both the 
emission and the excitation wavelengths. The refractive index of the medium and of the external 
is 1.4 at both the emission and the excitation wavelengths. An excitation spatial and temporal 
isotropic Dirac delta source of unitary strength is at (0, 0, l/fl s ). The absorption coefficient of 
the fluorophore at the excitation wavelength is fi a f = 0.002 mm -1 and 0 at the emission wave- 
length and its quantum efficiency is unitary. The same scattering properties as those for the 
background medium were present in the fluorophore. In Fig. 2 the TR reflectance fluorescence 
signal is obtained with the Homogeneous _Fluence .Fluorescence JSlab subroutine (that calculate 
Eq. (11) with the hybrid model for the slab) and with the partial current boundary condition. 
The hybrid model is calculated by using Paasschens RTE solution for the infinite medium [38]. 
The TR reflectance shown in Fig. 2 is the probability per unit time and per unit surface area 
that a photon emitted by the fluorophore exits from the medium at distance p and time t. In 
Fig. 2 we have assumed that fluorescence photons are emitted with a T =1 ns lifetime. The 
computation time for calculating Eq. (11) depends on the number of temporal intervals used 
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for calculating the time convolution integral. This parameter can be set to the right value in- 
side the Homogeneous -Fluence -Fluorescence Slab subroutine (Integer parameter denoted nd). 
By using 200 time windows (it is sufficient to have a negligible error in the calculated fluence 
due to the numerical evaluation of the integral): the computation time for calculating Eq. (11) 
for a single time is about 0.1 s with an Intel Centrino2 processor. In Fig. 3 a map of the CW 
reflectance fluorescence signal (Eq. (12)) is plotted that was obtained using the hybrid model 
(Homogeneous -Fluence -Fluorescence Slab subroutine) based on the IDE solution by Graaff 
and Rinzema and by using the same case as Fig. 2. The unit of the plotted CW reflectance is 
mm~ 2 , and represents the probability per unit surface area that a photon emitted by the fluo- 
rophore exits the from slab at a certain position on the z = 0 plane. The computation time of 
the whole map of Fig. 2 (mapping density of 1600 points) is about 0.01 s. 





X (mm) 



Fig. 3. Identical to Fig. 2 except for the Fluorescence CW reflectance. 

In Fig. 4 the TR reflectance fluorescence signal obtained with the Homoge- 
neous -Fluence -Fluorescence Slab subroutine and with the partial current boundary condition 
is compared with the reflectance simulated with a Monte Carlo code for a slab 2000 mm thick 
with the receiver at p =30 mm. For obtaining the Monte Carlo results we have modified a 
previously developed Monte Carlo code for describing the effects of absorbing and scattering 
defects [43]. The reduced scattering coefficient of the background medium is assumed equal 
to 0.5 mm -1 at both the emission and the excitation wavelengths. The fluorophore is assumed 
by a sphere centered at (15, 0, 15) mm. The same scattering properties as those for the back- 
ground medium were present in the fluorophore. In the MC simulations four different couples 
of values for the volume V' of the sphere and for the absorption jj. a f of the fluorophore at 
the excitation wavelength were considered (for each case we have V'jJ. a f = 0.01 mm ): V' = 
10 mm 3 , H a f = 0.001 mm" 1 ; V' = 100 mm 3 , \l af = 0.0001 mm" 1 ; V = 500 mm 3 , \i af = 
0.00002 mm" 1 , and V' = 1000 mm 3 , \l a} = 0.00001 mm -1 . The absorption coefficient of 
the fluorophore at emission wavelength is 0. The absorption coefficient of the background 
medium is 0.01 mm -1 in Fig. 4a, and 0 in Fig. 4b at both the emission and the excitation 
wavelengths. The refractive index of the medium and of the external is 1.4 at both the emission 
and the excitation wavelengths. In the comparison of Fig. 4 we have assumed that fluores- 
cence photons are promptly emitted. Thus in Eq. (11), and consequently inside the Homoge- 
neous -Fluence -Fluorescence Slab subroutine, a % =10~ 10 ps is set. The Monte Carlo results 
are obtained with an isotropic scattering function. The hybrid model is calculated by using in 
Eq. (11) and in the Homogeneous -Fluence -Fluorescence Slab subroutine the RTE solution for 
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the infinite medium. For the four different combinations of the volume V' and of the absorption 
jj, a f, the RTE hybrid model provides the same results since Eq. (11) depends only on the product 
V'lJ. a f that is constant. These results show that for the smallest volume V' the model is in pretty 
good agreement with the Monte Carlo results, while some deviations can be appreciated for 
the higher values of the volume V '. This fact confirms the expectations that the model provide 
a non sufficient description of photon migration for large volumes of the fluorophore. Similar 
results were also obtained by calculating the hybrid model for the slab with the DE solution for 
the infinite medium. 
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Fig. 4. Comparison between hybrid model and Monte Carlo results for the fluorescence TR 
reflectance signal received from a slab 2000 mm thick with the receiver at p = 30 mm. The 
reduced scattering coefficient of the background medium is assumed at the emission and 
the excitation wavelength equal to 0.5 mm -1 . The absorption coefficient of the background 
medium is 0.01 mm -1 in Fig. a), and 0 in Fig. b). In the MC simulations four different 
values of the volume and of the absorption of the fluorophore have been considered (for 
each case we have V'll a f = 0.01 mm 2 ): V' = 10 mm 3 , [l a f = 0.001 mm -1 ; V' = 100 mm 3 , 
j± af = 0.0001 mm" 1 ; V = 500 mm 3 , jx af = 0.00002 mm" 1 , and V' = 1000 mm 3 , \l af = 
0.00001 mm -1 . The fluorophore is assumed by a sphere centered at (15, 0, 15) mm. The 
refractive index of the medium and of the external is 1.4. The hybrid model is calculated 
with Eq. (11) and with the partial current boundary condition, and it has been implemented 
using the RTE solution for the infinite medium. 
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4. Summary and discussion 

We proposed a perturbative forward solver software for the time domain and the CW domain 
based on the Born approximation in order to simulate the effect of small localized fluorophores 
embedded inside a non-fluorescent biological tissue. The fluorescence emission is treated with a 
model that describes the migration of photons from the source to the fiuorophore and of emitted 
fluorescent photons from the fiuorophore to the detector, for all those geometries for which 
Green's functions can be calculated. The method has been implemented for the infinite medium 
and for the slab for which two subroutines (Homogeneous JFluence ^Fluorescence Infinite and 
Homogeneous -Fluence -Fluorescence Slab) written in FORTRAN language are provided with 
the linked file (Media 1). Quantities such as reflectance, transmittance, and fluence rate can 
be calculated using these subroutines and the software package distributed with a recently- 
published book [21]. Similar subroutines can be also obtained for other geometries such as the 
two-layer slab, the sphere, the parallelepiped, and the cylinder, by using the software distributed 
in [21] and by implementing for these geometries similar subroutines to those provided for 
the infinite medium and the slab. The software has been developed for the time and the CW 
domains, however it can also be extended to the frequency domain by following the guidelines 
provided in Sect. 2. The computation time of the said forward solver is significantly less than 
any elementary Monte Carlo program dedicated to the same kind of calculation. Typically, 
the computation time for calculating the TR reflectance (at a single value of time) emitted by a 
fluorescent inclusion is about 0.1s and the computation time for calculating the CW reflectance 
is less than 10~ 5 s. The accuracy of the proposed forward solver is excellent for small localized 
fluorescent sources, while it decreases when larger fluorescent inclusions are considered as it is 
shown in Fig. 4. The main use of forward solver software that we have proposed is for photon 
migration through biological tissues with embedded small localized fluorescent inclusions, such 
as those encountered in optical molecular imaging. 

The proposed forward solver is not meant to replace sophisticated numerical algorithms that 
are developed for solving partial differential equations in very general conditions. Nonetheless, 
we can apply this method to a certain number of geometries for which the Green's function 
can be expressed as a closed-form solution or it can be calculate numerically, and it has the 
advantage of being fast and precise for small localized fluorescent perturbations. In the case of 
larger fluorescent inclusions, it is possible to extend the validity of the theory to a wider range 
of fluorescent inclusions by implementing the perturbation theory beyond the limit of the Bom 
approximation [26-28,41,42]. In particular, the solutions obtained by Sassaroli et al. [28,41,42] 
for the higher order perturbation theory for CW, time and frequency domains have shown to be 
sufficiently accurate to describe correctly absorbing relative perturbations of intensities of up 
to 50-60%. For this reason, the forward solver that we propose could be further improved by 
using the perturbed solutions proposed in [28,41,42] for the Green's functions G and G e of 
Eqs. (11) and (12). 

In recent decades, several analytical and numerical solutions of the RTE have been reported. 
For instance, Machida et al. numerically obtained the Green's function for the RTE in the slab 
geometry for anisotropic scattering [39], and Liemert and Kienle obtained a CW analytical 
Green's function of the RTE for the infinite medium which is valid for anisotropic scattering 
[40]. We wish to emphasize that the availability of an analytical Green's function valid for 
anisotropic scattering makes it possible to further extend the perturbative approach presented 
in this paper. 
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